Velocity and diffusion coefficient of A + A <->• A reaction fronts in one dimension 
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We study front propagation in the reversible reaction-diffusion system A + A «-> A on a 1-d lattice. 
Extending the idea of leading particle in studying the motion of the front we write a master equation 
in the stochastically moving frame attached to this particle. This approach provides a systematic 
way to improve on estimates of front speed obtained earlier. We also find that the leading particle 
performs a correlated random walk and this correlation needs to be taken into account to get correct 
value of the front diffusion coefficient. 



I. INTRODUCTION 

Propagation of fronts into unstable states in reaction- 
diffusion systems has been an actively studied topic for 
a long time [1]. The dynamics of propagating fronts is 
of interest in various diverse physical situations [2]. Re- 
cently, there is renewed activity concerning connection 
between fronts in microscopic discrete stochastic models 
and macroscopic deterministic equations which are be- 
lieved to be their mean- field limits [3]. In this paper we 
study front propagation in a system of reacting and dif- 
fusing particles on an infinite 1-d lattice (Fig.l). Initially, 
left half of the lattice is filled with a certain density of 
particles with the right half being completely empty, re- 
sulting in a sharp boundary between the two halves. As 
the system evolves, the particles move to the right re- 
sulting in a propagating front. After an initial transient, 
the front moves with an asymptotic speed v. Further, 
due to the inherent stochastic nature of the dynamics, 
the ensemble averaged front profile undergoes diffusive 
broadening (la) with an associated diffusion coefficient 
Df- 




FIG. 1. Schematic picture of the development of front in 
the reaction-diffusion system, (a): The ensemble averaged 
front propagtes with a unique asymptotic speed v and spreads 
diffusively with a coefficient Df. (b) If the ensemble average 
is taken in a frame moving with the leading particle in each re- 
alization, the profile obtained is a time-invariant profile. The 
three profiles shown for different values of the parameters of 
the microscopic model discussed later in the text. 

There are a number of alternative definitons of the po- 
sition of the front in any given realization of the process 
but it is found that all of them lead to the same values of 
v and Df. The simplest choice, which is also the easiest 



to implement numerically, is to take the position of the 
rightmost particle at any time to be the location of the 
front [5-7]. Thus, the dynamics of the front is reduced, 
in an approximate way, to the dynamics of a single par- 
ticle - the leading particle. Treating the motion of this 
particle as a biased random walk, simple approximate 
expressions for the speed v and diffusion coefficient Df 
has been derived [7,6]. However, in these approaches it 
is hard to find a systematic way to improve estimates for 
v and Df . In the present paper we show that by writing 
the front dynamics in a frame moving with the leading 
particle one can get numerically better estimates of the 
front speed and diffusion coefficient. We also show that 
the motion of the leading particle is correlated in time 
and thus the front diffusion coefficient differs from that 
obtained from a simple random walk approximation. 



II. MODEL AND RESULTS FOR V AND D F 



We consider a 1-d lattice (-co < i < oo) in which 
each site can hold atmost one particle ('hard-core exclu- 
sion'). The particles (denoted by A) undergo three basic 
microscopic processes: (i) Birth/creation: a particle can 
generate a new one on a neighbouring empty site with 
rate e, (ii) Death/annihilation: one of the two neighbor- 
ing particles gets annhilated with rate W, and (iii) Dif- 
fusion: A particle diffuses to a neighbouring empty site 
with rate D (see Fig. 2). Initially, at t = 0, the left half 
of the lattice [i < 0) is filled with particles at a density 
p = ~p, where p = e/(e + W) is the density of the equilib- 
rium phase obtained if the process is allowed to occur in 
a finite system [8]. There are only two independent pa- 
rameters in the system as one of the three rates e,D,W 
can be scaled away by choosing the time scale appropri- 
ately. In all our simulations e = 0.25 and the two ratios 
D /e and W/ e are of the order of unity. As was noted in 
[7], this choice is fairly generic for the diffusion limited 
fluctuation dominated regime we would like to focus on. 
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FIG. 2. Basic microscopic moves of (a) diffusion, (b) birth 
and (c) death. When the right particle happens to be the 
leading particle (denoted by F) the processes can result in 
transitions between the two states Of and II as described 
later in the text. 

In the mean-field description, which is expected to 
be valid in the limit D — > oo and the evolution of the 
front is described by the F-KPPP equation [9] d t p — 
Dd 2 p + ep — (e + W)p 2 , where p is the coarse grained 
density of A particles. The front speed is then given by 
v = 2\[eD. 

For the case W = 0, which arises naturally in the con- 
text of turbulent flame front propagation [5], the exis- 
tence and uniqueness of an asymptotic front solution were 
established rigourously in [10] and it was shown that the 
'mean-field' limit is obtained as D — > oo. For finite val- 
ues of D, a 'two particle' representation was used in [6] to 
get an approximate value of the front speed which works 
quite well for D/e ~ 0(1). 

In [4], using the inter-particle distribution functions an 
exact solution was obtained for the special case W = D. 
In this case, it was shown that in each realization of the 
front evolution, the leading particle performs a biased 
random walk and the particles behind it are distributed 
exactly in equilibrium with density p = p. The speed 
and the diffusion coefficient were found to be v = e and 
Df — e + D respectively. 

If W ^ D, the method used in [4] does not work and an 
exact solution is no longer possible. Also, the two particle 
representation of [6] does not close if particle annihilation 
is introduced, i.e, W > 0. In [7], an approach based ex- 
plicitly on the motion of the leading particle as a biased 
random walker was introduced and the following approx- 
imate expressions for speed and diffusion coefficient were 
obtained. 

v = e- pi(W - D); 2D f = 2D + e - Pl (W - D), (1) 

where pi is the probability of occupation of the site just 
behind the leading particle. Expression (1) above may 
be written down by noting that the fronts moves right 
with rate P = e + D (i.e. whenever a particle is created 
to the right of the leading particle or the later makes a 
diffusive move to the right). The leading particle takes a 
negative step when it gets annihilated by the particle on 
the left (with rate Wpi) or if it makes a diffusive move 



to an empty left site (with rate D(l — pi)), where p\ is 
the occupation probability of the site immediately to the 
left of the leading particle. For W = D, these results 
reduce to those obtained exactly in [4]. In order to get 
more accurate values of v and Df one needs to find bet- 
ter estimates of p\ (the bulk value p\ = ~p was used in 
[7]). While (1) works quite well (with pi = p) for W close 
to D, there are significant deviations as one moves away 
from this special point. 

In order to get successively better estimates for p\, we 
look at the invariant profile of the front as observed from 
the leading particle. Starting from an ensemble of Af real- 
izations this invariant profile is obtained, after an initial 
transient time, by aligning the leading particle of each 
member of the ensemble (Fig. lb). From the definition 
of pi it then follows that out of the Af realizations p\M 
have a particle in the site to the left of the leading parti- 
cle and the rest (1 — p\)M have an empty site next to it. 
In the steady state (which is asymptotically reached after 
transients) it thus follows that there is a kinetic balance 
between the two types of realizations: those with 11 or 
01 as the occupancy of the rightmost pair of sites (where 
the second 1 denotes the leading particle). Thus, the 
two states 01 and 11 may be thought of as a truncated 
representation of the full lattice. Due to the microscopic 
moves there are transitions between these two 'states'. 
For example, in a realization in the 11 state the leading 
particle makes a diffusive move to the right, the state of 
the realization changes to 01 [11]. Considering all such 
transitions one can write a master equation for the prob- 
abilities pn and poi m this truncated state space [12] 



poi = (2D - D Pb + 2W) Pll - (2D Pb + 2e + ep b )p m , 
pn = (2D Pb + 2e + ep b )p i - {2D - Dp b + 2W) Pll . (2) 



In the steady state poi 
{Pu = Pi m steady state) 



= pn and one obtains 



Pi 



3e 2 + 2e(W + D) 



3e 2 + 2W 2 + 4eW + 3De + 2DW 



(3) 



where we have used p b — ~p, an approximation which 
becomes better as one includes more sites in the trun- 
cated representation. E.g. we have also computed p\ by 
keeping I = 3 sites (i.e. 4 states: 111, 011, 101, 001), 
/ = 4 sites (8 states) and the values obtained for p\ are 
plotted in Fig. 3 and the corresponding front speed (us- 
ing Eq.(l)) in Fig. 4. It is to be noted that expression 
(3) above reduces to p\ = p for W — D as it should. 
Further, in this approach it is also possible to obtain 
the spatial density correlation between site occupancies 
4>i2 = (nin 2 ) - P1P2/P 2 (Fig. 5). 



2 



0.014 
0.012 
0.01 h 
0.008 
0.006 
0.004 
0.002 


-0.002 
-0.004 
-0.006 



0.5 
0.4 
0.3 
0.2 
0.1 




— r~ 
1=1 
H=2 
l=3 
~3P 



i ■ * 



1=1 
l=2 
l=3 
3P 



» - 0.1 0.2 0.3 0.4*' 
■ ■ . 



0.05 



0.1 



0.15 



0.2 



0.25 



FIG. 3. Deviation of pi computed via various approxima- 
tions from that obtained from direct simulation as a function 
of D for W — 0.125 and W = (inset). Empty squares corre- 
spond to the 3-particle representation discussed in subsection 
A. 
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FIG. 4. Difference between front speed computed using pi 
from various approximations and that obtained from direct 
simulation as a function of D for W = 0.125 and W = 
(inset). Filled squares correspond to the 3-particle represen- 
tation discussed in subsection A. 



FIG. 5. The correlation 12 between occupancies of the 
pair of sites immediately following the leading particle as a 
function of D (W = 0.125). We see that the correlation is 
negative for D < W and positive for D > W. Inset: (f> 12 as a 
function of D for W = 0. 



A. Reduced 3-particle representation 

The drawback of the above approach is that it is not 
easy to write the transition matrix as the number of sites 
I is increased (number of states increases as 2 l ). In the 
following we try to find an analytically tractable esti- 
mate for p\ using an alternative reduced representation 
similar to Kerstein's approach [6]. Instead of keeping a 
fixed number of sites to denote a state, Kerstein chose 
the following infinite set of two particle states: {11, 101, 
1001, 10001, ••• } where the 'k'th state has k empty 
sites between the leading particle (denoted by the 1 on 
the right) and the next particle (denoted by the 1 on the 
left). If there is no annihilation process, i.e., W — 0, then 
these set of states is closed with respect to transitions be- 
tween the states. However, as was pointed out in [7], this 
breaks down as soon as W > 0. To see this let us con- 
sider the 11 state and the microscopic process in which 
one of the two particles gets annihilated. The resulting 
state depends upon the location of the third particle in 
the initial configuration of the lattice. I.e., one needs 
to go to a 3-particle representation. However the same 
problem occurs while considering annihilation in the 111 
state. Thus, to take care of this problem which arises due 
to the effective non-locality in transition rates by going 
to the moving frame attached to the leading particle one 
needs to make a reasonable truncation of the hierarchy. 
We extend the set of states by including the states {111, 
1011, 10011, • • • } where the rightmost 1 denotes the lead- 
ing particle. Let us denote by pk the probability of the 
two particle states 10^01 with k empty sites between 
k 

the leading particle and the next and by qk the proba- 
bility of the three particle state 10. ..011 with k empty 
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sites between the second and the third particles. The 
transition rates within these states can be written down 
following the procedure as in the case of fixed number of 
sites representation discussed above and one obtains the 
following rate equations for the pit's. 

p k = (2D -Dp + Wp)p k -i + (2D + e)p k+ i 

+ (q k -i + Qk)W - (AD -Dp + 3e + Wp)p k ; k>2, 
pi = (2D - Dp) Po + (2D + e)p 2 + (2q + qi )W 

-(4D- Dp + 3e + Wp) Pl ; 
po = (2D + t)pi + 2e(l - po) - (2D -Dp- 2W)p a , (4) 

where we made the approximation that p is the proba- 
bility that the site next to the last particle in the two or 
three particle states is occupied independent of its dis- 
tance from the leading particle. We do not write the 
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equations for q^s as it can be shown that they can be 
eliminated in the steady state. To obtain the steady 
state solution pk = 0, following Kerstein, we make the 
ansatz pk = po(l — Po) k and write p = Ap a — Bp\ to take 
care of the ignored correlations in the ansatz in a phe- 
nomenological way. The parameters A and B are fixed 
in the following way. We first note that for D = W, both 
the quantities p and po equal the bulk density namely 
p = Po = p = e /(e + W) implying A = l + eB/(e + W) 
and thus 



P=(l + 



eB 
e + W 



)Po - Bpl 



(5) 



This form also ensures that W — and pa=l imply p=l. 
Using p from Eq. (5) in Eqs. (4) (in steady state p\ = 0) 
we get the following cubic equation for po, 

DB(e + W)pl + (e 2 + eW + eD + WD - DBe)p 2 Q 
+ (e 2 + ieW + 2W 2 )p Q - 2e 2 - 2eW = (6) 

In order to fix B we note that for large D we expect 
the front to approach the mean-field limit with the speed 
given by vo = 2\feD. Since, in terms of po the velocity 
is given by (from eq. (1)) v = e — (W — D)po <~ Dpo 
for D >> W,e, this implies, in this limit, po = 2^/2e/D. 
Substituting this expression for p in Eq. (6) we obtain 
B = 3(e + W)/Ae. Using this expression for B in Eq. (6) 
one obtains po = pi implicitly through, 



D 



4e(2e - ep 2 - ep - 2 Po W) 
(e + 3ep + 3Wp )pl 



(7) 



The result is shown in Fig. 3 (square symbols marked 
as 2P). We note that, although several ad-hoc assump- 
tions were made in arriving at eq. (7), including us- 
ing results that are strictly valid in the mean-field limit 
D — > oo, the agreement with direct numerical results is 
remarkable even for the finite values of D considered. 



B. Diffusion coefficient: Df 

In Fig. 6, we plot the front diffusion coefficient Df as a 
function of D (W = 0, e = 0.05). The lower curve is from 
direct simulation and the upper one is that obtained from 
Eq. (1) by using the most accurate estimate of p\. We 
see marked deviation of the values obtained from the an- 
alytic expression which implies that the simple (uncorre- 
cted) random walk picture of the leading particle in not 
quite correct. The motion of a correlated biased random 
walker is described by x(t + 1) — x{t) = v + rj{t) where 
the noise term ij is temporally correlated: (77(f)) = and 
(v(t)v(t')) = C(f-f') ^ 0. The mean speed of the walker 
is v and the asymptotic diffusion coefficient is given by 



Indeed, for the front under study we find that there is 
long range correlation between the successive steps of 
the leading particle (Fig. 6, inset). This correlation is 
non-positive for all parameters (both D > W as well as 
D < W) and vanishes for the special case of D = W. 
Once this correlation is taken into account the diffusion 
coefficient matches reasonably with that obtained from 
direct simulations for the range of D studied [13]. Pre- 
liminary fits indicate that the correlation function has 
the functional form C(f) = A t~ a exp (— f/r). 
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FIG. 6. Diffusion coefficient of the front for W = as a 
function of D. The bottom data (pluses) are the direct simu- 
lation values. The top (crosses) are the values obtained from 
Eq. (1) while the middle one (stars) represent the correlation 
corrected Df. Inset: the correlation function C(f) for differ- 
ent values of D. 



III. CONCLUSION 

We have illustrated the usefulness of the leading par- 
ticle picture in describing the propagation of fronts in 
the A + A <-» A reaction-diffusion process in the diffusion 
controlled limit in one dimension. By writing the mas- 
ter equation in the moving frame attached to the leading 
particle we are able to obtain better numerical estimates 
for the density of the site behind the leading particle and 
thus the front speed. In addition, this approach, in prin- 
ciple allows one to compute the spatial density profile and 
density-density correlations away from the special point 
D = W. 

Our numerical results show that the motion of the lead- 
ing particle is correlated in time and this needs to be 
taken into account in order to get the correct diffusion 
coefficient. It is seen that this correction increases with 
increasing D (the microscopic particle diffusion constant) 
and thus it might play an important role in determining 
how the mean- field limit is achieved as D — > 00. 



D / = ^(77(0)77(f))=£c(r). 



(8) 
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